3 Stochastic Differential Equations – Advanced Mathematical Biology IV (MATH4411)

In this chapter we introduce a new class of continuous-time continuous-state stochastic processes that are commonly used for modelling in mathematical biology. Besides being of independent interest, this class of processes will offer more analytic tractability than the chemical reaction processes of the first two chapters of the course. As we have seen in examples and problem classes, analytic calculations for chemical reactions processes with many species and reactions, particularly higher-order reactions, quickly become impossible with just pen and paper! One advantage of these processes is that we can use them to approximate chemical reactions processes and thereby study the complex noise-induced phenomena of Chapter 2 in more detail. The new stochastic processes we introduce will have dynamics described by so-called stochastic differential equations (SDEs), which we can think of as generalizations of ordinary differential equations with random noise added. Instead of a fully rigorous treatment of SDEs, which is an entire course in its own right, we will proceed by discretising time and introducing a “computational definition” of SDEs that gives us the appropriate intuition, enables basic calculations and allows us to perform numerical simulations that produce the correct dynamics.

3.1 Introduction to SDEs

Consider the ordinary differential equation

\[ \frac{d}{dt}x(t) = f(x(t)), \quad t >0, \quad x(0)= x_0. \tag{3.1}\]

We can also write Equation 3.1 in the form

\[ dx(t) = f(x(t))\,dt, \]

where we interpret \(dx\) as an increment in \(x\) and \(dt\) as a (small) increment in time. In this context, we mean \(dx(t) \approx x(t+dt) - x(t)\) and this naturally leads to the discretisation of the ODE Equation 3.1,

\[ x(t + \Delta t) = x(t) + f(x(t))\Delta t, \tag{3.2}\]

where we choose \(\Delta t>0\) as the discretisation parameter. Hence, given the initial condition, we can use this discretisation to numerically approximate the solution to Equation 3.1 by stepping forward in time in steps of \(\Delta t\) (the forward Euler scheme). If \(f\) is sufficiently smooth and we choose \(\Delta t\) sufficiently small, we obtain a good approximation and if we let \(\Delta t \downarrow 0\), we recover back the differential equation Equation 3.1 exactly (formally, the scheme is first order convergent in \(\Delta t\)).

We could think of Equation 3.2 as a discrete time stochastic process. However, as stated, it is totally deterministic as there is no source of randomness… but we could add some! For example, we could write

\[ x(t + \Delta t) = x(t) + f(x(t))\Delta t + G(t), \]

where \(G(t)\) is some stochastic process. For reasons that will become clear presently (and some others that won’t), we need to insist on more structure than this; we instead choose to extend Equation 3.2 to processes of the form:

\[ X(t + \Delta t) = X(t) + f(X(t))\Delta t + g(X(t)) \sqrt{\Delta t}\,\xi, \tag{3.3}\]

where \(\xi\) is a normally distributed random variable with mean \(0\) and variance \(1\). Every time we take a step \(\Delta t\) forward in time, we draw a new realisation of \(\xi\), so we really have a discrete collection of standard normal random variables generating the “noise” in the process Equation 3.3. Fortunately, all modern programming languages can generate normal random variables and hence this makes processes such as Equation 3.3 straightforward to simulate on a computer. Note that if \(g\equiv 0\), then Equation 3.3 reduces exactly to Equation 3.2. If \(g\) is in some sense small, we can also expect that the dynamics of \(X\) will often closely approximate those of the solution to Equation 3.2 (although there will be notable exceptions to this intuition!). Finally, when we write Equation 3.3, we typically restrict \(t\) to the discrete set \(\{ \Delta t, 2\Delta t, \dots \}\) and we can interpolate linearly between these points if we want to know the value of \(X(t)\) between time points.

We will refer to equation Equation 3.3 as our computational definition of stochastic differential equations (SDEs).

The SDE Equation 3.3 is more often expressed in the form:

\[ dX(t) = f(X(t))\,dt + g(X(t)) \, dW(t), \tag{3.4}\]

where \(W(t)\) is a process whose increments, i.e., \(dW(t)\), are normally distributed. In fact, both Equation 3.3 and Equation 3.4 are formal representations of the process \(X(t)\). Rigorously, the equation describing the dynamics of \(X(t)\) is actually the integral equation

\[ X(t) = X(0) + \int_0^t f(X(s))\,ds + \int_0^t g(X(s)) \, dW(s), \tag{3.5}\]

which can be formally obtained by integrating Equation 3.4. We can think of Equation 3.5 as equivalent to writing the ODE Equation 3.1 as \(x(t) = x(0) + \int_0^t f(x(s))\,ds\) by simply integrating across with respect to \(t\). The final term on the right-hand side of Equation 3.5 is technically a “stochastic integral” and defining such an object rigorously requires developing a whole new theory of integration (which is well beyond the scope of this course). In numerical analysis terms, equation Equation 3.3 is referred to as the Euler-Maruyama scheme for discretising Equation 3.5.

You don’t need to worry about the details of Equation 3.5 but it is good to be aware of as you may see the same SDE written in different ways in other contexts!

We can immediately write down an SSA for Equation 3.3:

At time \(t = 0\), set \(X(0)\), then:

Generate \(\xi \sim N(0,1)\).
Set \[ X(t + \Delta t) = X(t) + f(X(t))\Delta t + g(X(t)) \sqrt{\Delta t}\,\xi, \] update \(t = t+\Delta t\), and go back to step 1.

Example 3.1

Choose \(X(0)=0\), \(f \equiv 0\) and \(g \equiv 1\). Then Equation 3.3 reads:

\[ X(t + \Delta t) = X(t)+ \sqrt{\Delta t}\,\xi. \tag{3.6}\]

Suppose \(M(t) = \mathbb{E}[X(t)]\) and \(V(t) = \text{Var}[X(t)]\). We can then calculate as follows using our computational definition:

\[ \begin{aligned} M(t+\Delta t) &= \mathbb{E}[X(t)] + \mathbb{E}\left[\sqrt{\Delta t}\,\xi \right] = M(t). \end{aligned} \]

Since \(X(0)= 0\), \(M(0) = 0\) and thus \(M(t) = 0\) for all \(t \geq 0\). Similarly, we have

\[ \begin{aligned} V(t+ \Delta t) &= \mathbb{E}\left[ X(t+\Delta t)^2 \right] - M(t + \Delta t)^2\\ &= \mathbb{E}\left[ X(t) + \sqrt{\Delta t}\, \xi)^2 \right] \\ &= \mathbb{E}[X(t)^2] + 2 \sqrt{\Delta t}\mathbb{E}[X(t)] \mathbb{E}[\xi] + \Delta t \mathbb{E}[\xi^2]\\ &= \mathbb{E}[X(t)^2] + \Delta t\\ &= V(t) + \Delta t. \end{aligned} \]

Since \(V(0) = 0\), \(V(t) = t\) for all \(t \geq 0\). In other words, the variance of the process Equation 3.6 scales linearly with time. This is a crucial aspect of our choice in adding the noise term when we wrote down Equation 3.3. Moreover, this choice has ensured that both the mean and variance of the process are independent of the time step, suggesting that this was somehow the “right” scaling to choose for the noise term. In fact, we can show that all moments of the process Equation 3.6, i.e., \(\mathbb{E}[X(t)^k]\), are independent of step size \(\Delta t\).

Sample paths of the process computed via the SSA above are shown in Figure 3.1 (left panel). We observe that these paths are quite centred around zero with larger fluctuations away from zero as time increases, as expected given that variance scales linearly with time for the process Equation 3.6.

Example 3.2

We can consider higher dimensional examples of SDEs, such as the process

\[ \begin{aligned} X(t+ dt) &= X(t) + \sqrt{2D}\, dW_x, \\ Y(t+ dt) &= Y(t) + \sqrt{2D}\, dW_y, \\ Z(t+ dt) &= Z(t) + \sqrt{2D}\, dW_z, \end{aligned} \tag{3.7}\]

where \(D > 0\) is a parameter and \(W_x\), \(W_y\) and \(W_z\) are three independent sources of noise. For simulations, as usual, we replace \(dW_x\) by \(\sqrt{\Delta t}\, \xi\) where \(\xi \sim N(0,1)\) at each time step. The process Equation 3.7 may be interpreted as a stochastic model for a diffusing particle in 3 spatial dimensions, with \(D>0\) denoting the diffusion coefficient. The right-hand panel of Figure 3.1 below shows some sample paths of Equation 3.7 using the SSA outlined above with \((X(0),Y(0),Z(0)) = (0,0,0)\) for all simulations. The end point of each simulation is marked with a black dot and as we follow each path, at least qualitatively, we observe what we expect to from a particle diffusing freely in an aqueous medium.

Figure 3.1: Left: Sample paths of the SDE Equation 3.6. Right: Sample paths of the process Equation 3.7 projected onto the \(X\)-\(Y\) plane. Parameters: \(D = 0.0001\).

Exercise 3.1

Open the course Github page and try playing with the MATLAB scripts

CH3_example_1.m
CH3_example_2.m

to see how the dynamics vary as you increase the discretisation parameter and the diffusion coefficient.

Example 3.3

Choose \(f \equiv 1\) and \(g \equiv 1\) so that Equation 3.3 reads

\[ X(t + \Delta t) = X(t) + \Delta t + \sqrt{\Delta t}\,\xi = X(t) + dW(t), \quad X(0)=0. \tag{3.8}\]

In this case, we can compute as in Example 1 to show that \[ M(t + \Delta t) = M(t) + \Delta t \implies M(t) = t, \] using that \(M(0) = 0\) once more.

Figure 3.2 (left) shows some sample paths of the process Equation 3.8. The black dashed line denotes the function \(M(t) = t\) and we see that all of the sample paths roughly follow this trend line. This trend line is solely set by the function \(f \equiv 1\) which is typically referred to as the “drift” coefficient of the SDE. The function \(g\) is called the “diffusion” coefficient and controls the noise level in the process.

Example 3.4

For our final example, choose \(f(x) = -k_1 x^3 + k_2 x^2 -k_3 x + k_4\) and \(g(x) = k_5\).

Hence we have the SDE \[ X(t + \Delta t) = X(t) + \left( -k_1 X(t)^3 + k_2 X(t)^2 -k_3 X(t) + k_4 \right)\Delta t + k_5\sqrt{\Delta t}\,\xi. \tag{3.9}\]

If we set \(k_5 = 0\), Equation 3.9 is exactly the ODE we studied in Section 2.1 when we used the law of mass action to derive a deterministic model for the bistable chemical reaction process Equation 2.1. The mass action model for the process Equation 2.1 was given by

\[ \frac{d}{dt} a(t) = -k_1 a(t)^3 + k_2 a(t)^2 -k_3 a(t) + k_4, \quad t \geq 0. \tag{3.10}\]

If we choose reaction rate parameters

\[ k_1 = 0.001, \quad k_2 = 0.75, \quad k_3=165, \quad k_4=10000, \] then the ODE Equation 3.10 has stable steady states \(\bar{A}_1 = 100\) and \(\bar{A}_3 = 400\), and an intermediate unstable steady state \(\bar{A}_2 = 250\). We saw that the solution to the ODE is not a good model or approximation of the bistable chemical reaction because it cannot switch between the stable fixed points, as the stochastic reaction process can.

In Figure 3.2 we show simulations of the process Equation 3.9 with \(k_5 = 200\) and observe that the solution to the SDE can also switch between the stable fixed points of the deterministic model. This illustrates that SDEs offer the potential of better qualitatively matching the dynamics of stochastic reaction processes. Moreover, we can tune the noise level in the SDE to understand how much noise is required to induce switching dynamics and, as we will see presently, we can use the analytic tractability of the SDEs to understand the mean switching time in chemical reaction systems.

Figure 3.2: Left: Sample paths of the SDE Equation 3.8. Right: Sample path of the process Equation 3.9 and solutions to the corresponding ODE with zero noise.

Exercise 3.2

Open the course Github page and use the MATLAB script

CH3_example_4.m

to investigate how the switching dynamics are impacted by the noise level in Example 3.4 by adjusting the parameter \(k_5\).

3.2 The Fokker-Planck Equation

3.2.1 Overview

When we studied chemical reaction processes, we derived the chemical master equations, which can be solved to yield the probability mass function of the process at each point in time. Similarly, if the process \(\{X(t) \, :\, t \geq 0\}\) solves Equation 3.3 and we denote its probability density function by \(p(x,t)\), we can derive an evolution equation for \(p(x,t)\).

Since \(p(x,t)\) is a PDF, we must have

\[ \int_\mathbb{R} p(x,t)\, dx = 1, \]

and we have the usual intuitive interpretation of \(p(x,t)\) as a proxy for the probability that \(X(t)\) is near \(x\) at time \(t\), i.e.

\[ \mathbb{P}[X(t) \in [x,+dx]] \approx p(x,t) \, dx. \]

We will show that \(p(x,t)\) solves the linear second-order PDE:

\[ \boxed{\frac{\partial }{\partial t} p(x,t) = \frac{\partial^2}{\partial x^2}\left( \frac{1}{2}g^2(x) p(x,t) \right) - \frac{\partial }{\partial x}\left( f(x)p(x,t) \right), \quad t > 0.} \tag{3.11}\]

Equation 3.11 is typically called the Fokker-Planck equation across most fields, although it can also be referred to as the Kolmogorov forward equation or the Smoluchowski equation.

As an example, consider the case when \(f \equiv 0\) and \(g \equiv \sqrt{2 D}\) for some \(D>0\). Then Equation 3.11 reads

\[ \frac{\partial }{\partial t} p(x,t) = D \frac{\partial^2}{\partial x^2} p(x,t). \]

In other words, with no drift coefficient and constant diffusion, the Fokker-Planck equation is just the heat equation. Therefore the solution in this case is

\[ p(x,t) = \frac{1}{\sqrt{2 \pi t D}} e^{-x^2/ 2D t}, \quad t \geq 0, \]

if we suppose that the initial condition for the process is \(X(0) = 0\); this initial condition corresponds to solving the PDE with the initial distribution being a Dirac delta function (i.e., all mass is at \(x = 0\) for \(p(x,0)\)). Figure 3.3 (left) shows the solution to the Fokker-Planck equation for this example, along with the estimated PDF from simulations at time \(t = 1\). As expected, agreement between the formula and the SSA approach is very close.

When the drift and diffusion coefficients, \(f\) and \(g\), are independent of time, we can characterize the long-term behaviour of the SDE Equation 3.3 in terms of the stationary distribution of the process, \(p_s(x)\), i.e.

\[ p_s(x) = \lim_{t \to \infty} p(x,t). \]

It follows that \(p_s(x)\) obeys the steady-state version of the PDE Equation 3.11:

\[ 0 = \frac{1}{2}\frac{d^2}{dx^2}\left( g^2(x) p_s(x) \right) - \frac{d}{dx}(f(x) p_s(x)). \tag{3.12}\]

In problem classes, we showed that the solution to Equation 3.12 is formally given by

\[ \boxed{p_s(x) = \frac{C}{g^2(x)} \exp\left( \int_0^x \frac{2f(y)}{g^2(y)}\,dy \right),} \tag{3.13}\]

where the normalisation constant \(C>0\) is given by

\[ C = \left(\int_{\mathbb{R}}\frac{1}{g^2(x)} \exp\left( \int_0^x 2f(y)/g^2(y)\,dy \right)\,dx\right)^{-1}. \]

Figure 3.3: Left: Comparison of the Fokker-Planck solution and estimated probability density with \(f \equiv 0\) and \(g \equiv 1\) at time \(t = 1\). Right: Estimated stationary distribution from simulations and stationary distribution obtained via the Fokker-Planck equation for the bistable process Equation 3.9. [Codes: `CH3_FP1.m`, `CH3_FP2.m`]

We saw in our first example above that when drift and diffusion are constant, we obtain a Normal distribution for the process at all times. However, for different choices, both the finite-time and stationary distributions generated can be quite far from Normal. For example, Figure 3.3 (right) shows a comparison of the stationary distribution for the bistable process Equation 3.9 from the last section obtained via both simulations and solving the Fokker-Planck equation using the formula Equation 3.13. Both approaches produce a bimodal distribution in this case, reflecting the stochastic switching between the two stable states at \(x = 100\) and \(x = 400\) (see Figure 3.2).

3.2.2 Derivation

To derive the Fokker-Planck equation Equation 3.11, we begin by defining the conditional probability density \(p(x,t \,|\, y,s)\), which is approximately the probability that \(X(t) \in [x,x+dx]\) given that \(X(s) = y\) for \(s<t\).

Now consider the value of \(X(t + \Delta t)\) for some \(\Delta t > 0\). We can divide the interval \([s,t+ \Delta t]\) into the disjoint intervals \([s,t]\) and \((t,t+\Delta t]\). If we want to arrive at \(X(t + \Delta t) = z\), we can think of all possible paths through an intermediate point \(X(t) = x\). In other words, we may write

\[ \boxed{p(z,t+\Delta t \,|\, y,s) = \int_{\mathbb{R}} p(z,t+\Delta t \,|\, x,t) \,p(x,t \,|\, y,s)\,dx, \quad s < t.} \tag{3.14}\]

Equation 3.14 is called the Chapman-Kolmogorov equation and is valid for all \(\Delta t > 0\), even though our goal here will be to later let \(\Delta t \downarrow 0\) in order to obtain an evolution equation for \(p(x,t \,|\, y,s)\).

Next, to alleviate potential regularity issues, we multiply Equation 3.14 across by a smooth test function \(\phi(z)\) and integrate over \(z\) to obtain: \[ \int_{\mathbb{R}} p(z,t+\Delta t \,|\, y,s) \phi(z)\,dz = \int_{\mathbb{R}} \left[\int_{\mathbb{R}} p(z,t+\Delta t \,|\, x,t)\phi(z)dz \right] \,p(x,t \,|\, y,s)\,dx, \] which in the long run will be less confusing to rewrite as \[ \int_{\mathbb{R}} p(x,t+\Delta t \,|\, y,s) \phi(x)\,dx = \int_{\mathbb{R}} \left[\int_{\mathbb{R}} p(z,t+\Delta t \,|\, x,t)\phi(z)dz \right] \,p(x,t \,|\, y,s)\,dx, \tag{3.15}\]

where we have simply changed the integration variable on the left-hand side. Now Taylor expand the (smooth test function) \(\phi\) about the point \(z = x\) on the right-hand side of Equation 3.15:

\[ \begin{aligned} &\int_{\mathbb{R}} \left[\int_{\mathbb{R}} p(z,t+\Delta t \,|\, x,t)\phi(z)dz \right] \,p(x,t \,|\, y,s)\,dx = \\ &\quad \int_{\mathbb{R}}\int_{\mathbb{R}} p(z,t+\Delta t \,|\, x,t) \left\{ \phi(x) + \phi '(x)(z-x) + \phi''(x) \tfrac{(z-x)^2}{2} + o((z-x)^2) \right\} p(x,t \,|\, y,s)\,dx. \end{aligned} \tag{3.16}\] Now we need to tackle some of the integrals on the right-hand side. Happily, we can rewrite these as (conditional) expectations and compute them using the computational definition of the underlying SDE, i.e., Equation 3.3. Firstly, we have \[ \begin{aligned} \int_{\mathbb{R}} (z-x) p(z,t+\Delta t \,|\, x,t) \, dz &= \mathbb{E}\left[ X(t + \Delta t) - x \, | \, X(t) = x \right] \\ &= \mathbb{E}\left[ f(x)\Delta t + g(x) \sqrt{\Delta t}\, \xi \right] = f(x)\Delta t. \end{aligned} \tag{3.17}\]

Similarly, we can compute the second integral as:

\[ \begin{aligned} \int_{\mathbb{R}} (z-x)^2 p(z,t+\Delta t \,|\, x,t) \, dz &= \mathbb{E}\left[ \left(X(t + \Delta t) - x\right)^2 \, | \, X(t) = x \right] \\ &= \mathbb{E}\left[ \left( f(x)\Delta t + g(x) \sqrt{\Delta t}\, \xi \right)^2 \right] \\ &= \mathbb{E} \left[ f^2(x) \Delta t^2 + 2 f(x) g(x) (\Delta t)^{3/2} \xi + \xi^2 g^2(x) \Delta t \right] \\ &= g^2(x)\Delta t + o(\Delta t^2). \end{aligned} \tag{3.18}\]

We have truncated at order \(\Delta t\) since we intend to take \(\Delta t \downarrow 0\) later and these higher order terms will vanish; this is also why we truncated our expansion of \(\phi\) at second order. Returning to Equation 3.15 with these estimates and plugging them in yields: \[ \int_{\mathbb{R}} p(x,t+\Delta t \,|\, y,s) \phi(x)\,dx = \int_{\mathbb{R}} \left\{ \phi(x) + \phi '(x)f(x)\Delta t + \phi''(x) \frac{g(x)^2}{2} \Delta t \right\} p(x,t \,|\, y,s)\,dx + o(\Delta t^2). \] With some rearrangement and division by \(\Delta t\) on both sides, we thus have \[ \begin{aligned} \int_{\mathbb{R}} \frac{p(x,t+\Delta t \,|\, y,s) - p(x,t \,|\, y,s)}{\Delta t} \phi(x)\,dx = \int_{\mathbb{R}} \left\{\phi '(x)f(x) + \phi''(x) \frac{g(x)^2}{2} \right\} p(x,t \,|\, y,s)\,dx + o(\Delta t). \end{aligned} \tag{3.19}\] The integrals on the right-hand side of Equation 3.19 can be simplified using integration by parts. Under mild regularity assumptions on \(f\) and \(g\), and noting that \(\lim_{|x| \to \infty}p(x,t \, | \, y,s) = 0\), we obtain

\[ \begin{aligned} \int_{\mathbb{R}} \phi '(x)f(x) p(x,t \,|\, y,s)\,dx &= - \int_{\mathbb{R}}\phi(x) \frac{\partial }{\partial x}\left( f(x) p(x,t \,|\, y,s) \right)dx \\ \int_{\mathbb{R}} \phi''(x) \frac{g(x)^2}{2} \,p(x,t \,|\, y,s)\,dx &= \int_{\mathbb{R}} \phi(x)\frac{\partial^2 }{\partial x^2}\left( \frac{g^2(x)}{2} p(x,t \,|\, y,s) \right)dx. \end{aligned} \]

Using these identities in Equation 3.19 and letting \(\Delta t \downarrow 0\) thus yields

\[ \int_{\mathbb{R}} \phi(x) \frac{\partial }{\partial t} \, p(x,t \,|\, y,s) \,dx = \int_{\mathbb{R}} \phi(x)\left\{- \frac{\partial }{\partial x}\left( f(x) p(x,t \,|\, y,s) \right) + \frac{\partial^2 }{\partial x^2}\left( \frac{g^2(x)}{2} p(x,t \,|\, y,s) \right) \right\}dx, \]

or equivalently,

\[ \int_{\mathbb{R}} \phi(x)\left\{\frac{\partial }{\partial t} \, p(x,t \,|\, y,s) - \frac{\partial }{\partial x}\left( f(x) p(x,t \,|\, y,s) \right) + \frac{\partial^2 }{\partial x^2}\left( \frac{g^2(x)}{2} p(x,t \,|\, y,s) \right) \right\}dx = 0. \]

The only way this equality can be satisfied for an arbitrary test function \(\phi\) is if the integrand is identically zero, i.e.

\[ \frac{\partial }{\partial t} \, p(x,t \,|\, y,s) = - \frac{\partial }{\partial x}\left( f(x) p(x,t \,|\, y,s) \right) + \frac{\partial^2 }{\partial x^2}\left( \frac{g^2(x)}{2} p(x,t \,|\, y,s) \right). \tag{3.20}\]

The PDE Equation 3.20 is valid for any time \(s<t\) and any \(y\) but for the initial value problem for the SDE Equation 3.3, we simply want the case \(s = 0\) and \(y = x_0\). Thus we can let \(p(x,t) = p(x,t \,|\, x_0,0)\) and specialize Equation 3.20 to obtain \[ \boxed{ \frac{\partial }{\partial t} \, p(x,t) = - \frac{\partial }{\partial x}\left( f(x) p(x,t) \right) + \frac{\partial^2 }{\partial x^2}\left( \frac{g^2(x)}{2} p(x,t ) \right),} \] which is exactly the claimed Fokker-Planck equation given by Equation 3.11.

3.2.3 Boundary conditions

In the previous section we derived the Fokker-Planck equation under the assumption that the underlying process \(X\) took values on all of \(\mathbb{R}\). In effect, there is no boundary and hence there is no need for a boundary condition to solve Equation 3.11 in this scenario.

In many practical situations we want to restrict the set of values that the process \(X\) can attain and this has natural implications for the Fokker-Planck equations of the process. For example, if we are modelling the positions of particles or cells in a Petri dish, then the domain of the cell positions is probably best modelled as a compact subset of \(\mathbb{R}^2\). We must then decide what happens to cells when they come into contact with the boundary of the dish.

Suppose that a particle moves in 1 dimension and its position must remain above the level \(x = 0\). We can impose a “reflecting” boundary condition at \(x = 0\) to achieve this, i.e., the particle hits the boundary and is pushed back into the admissible domain \((0,\infty)\). If we are using the recurrence relation Equation 3.3 to simulate the process, we can modify our previous SSA for SDEs as follows to account for the addition of a reflecting boundary:

At time \(t = 0\), set \(X(0)\), then:

Generate \(\xi \sim N(0,1)\).
Compute the possible value of \(X(t+\Delta t)\): \[ X(t + \Delta t) = X(t) + f(X(t))\Delta t + g(X(t)) \sqrt{\Delta t}\,\xi, \]
If \(X(t+\Delta t) \geq 0\), then we accept that value. If \(X(t+ \Delta t) < 0\), then we instead take \[ X(t + \Delta t) = - X(t) - f(X(t))\Delta t - g(X(t)) \sqrt{\Delta t}\,\xi. \] Finally, set \(t = t+\Delta t\) and go back to step 1.

In the “reflective step” in step 3 above, we must have \(X(t)> 0\). Thus if \(X(t+\Delta t) < 0\) in step 2, we must have had

\[ f(X(t))\Delta t + g(X(t)) \sqrt{\Delta t}\,\xi < 0. \]

In the reflective step, we first reflect the process through zero by replacing \(X(t)\) by \(-X(t)\). Then we add the positive increment \(- f(X(t))\Delta t - g(X(t)) \sqrt{\Delta t}\,\xi\) to make the process positive once more.

As an example, consider the SDE Equation 3.3 with \(f(x) = -x\) and \(g \equiv 1\) with initial condition \(X(0) = 0\). We can simulate a path of this process using our SSA with and without the reflecting boundary condition to see the impact of the boundary on the process. To make the comparison even more direct, we will use the same set of normal random variables (\(\xi\)’s) for the noise terms. In the left-hand panel of Figure 3.4 we see how the paths of the reflected and unreflected processes diverge as soon as they hit zero for the first time. The reflected process (blue) hits zero numerous times but is blocked from crossing into \((-\infty,0)\), while the red process is unrestricted. We observe that both processes have the same dynamics away from \(x = 0\) since they share the same underlying noise process and thus their increments are mostly the same (except for reflection events). The right-hand panel of Figure 3.4 shows the estimated PDF of the reflected and unreflected processes at time \(t =10\). As expected, the PDF of the reflected process is only supported on \([0,\infty)\).

Figure 3.4: Left: Samples path of \(X(t)\) and \(X_R(t)\) with \(f(x) = -x\) and \(g \equiv 1\) with identical initial conditions and noise processes (i.e. the same random variables \(\xi\)). Right: Estimated PDFs for \(X\) and \(X_R\) at time \(t = 10\). [Code: `CH3_FP3_reflecting.m`]

Figure 3.4 shows how significantly boundary conditions can change the solution to the Fokker-Planck equation, this immediately begs the question: How does the reflecting boundary condition impact the derivation of the corresponding Fokker-Planck equation?

Let \(p_R(x,t\,|\, y,s)\) denote the conditional probability density function of a process \(X_R\) solving Equation 3.3 subject to the reflecting boundary condition introduced above. The process \(X_R\) can no longer take on negative values with positive probability so \(p_R(x,t \, | \, y,s) = 0\) if \(x<0\) or \(y<0\). The Chapman-Kolmogorov equation Equation 3.14 then becomes

\[ p_R(z,t+\Delta t \,|\, y,s) = \int_{[0,\infty)} p_R(z,t+\Delta t \,|\, x,t) \,p_R(x,t \,|\, y,s)\,dx, \quad s < t. \tag{3.21}\]

As before, we can multiply by a smooth test function \(\phi\) and integrate to obtain:

\[ \int_{[0,\infty)} p_R(x,t+\Delta t \,|\, y,s) \phi(x)\,dx = \int_{[0,\infty)} \left[\int_{[0,\infty)} p_R(z,t+\Delta t \,|\, x,t)\phi(z)dz \right] \,p_R(x,t \,|\, y,s)\,dx. \tag{3.22}\]

Next consider the unreflected process \(X\), whose conditional density is denoted \(p(x,t \, | \, y,s)\), and whose Fokker-Planck equation we derived in the previous section. We assume that \(X\) solves the same SDE as \(X_R\). If \(x>0\), \(y>0\), then a path of \(X\) that goes from \(X(s) = y\) to \(X(t) = -x\), has the same probability as a path of the reflected process that goes from \(X_R(s) = y\) to \(X_R(t) = x\) due to the symmetry of the reflection step. If our test function \(\phi\) is chosen to be an even function, we then have

\[ \begin{aligned} \int_{[0,\infty)} p_R(z,t+\Delta t \,|\, x,t)\phi(z)dz &= \int_{[0,\infty)} p(z,t+\Delta t \,|\, x,t)\phi(z)dz \\ &\qquad + \int_{(-\infty,0]} p(z,t+\Delta t \,|\, x,t)\phi(-z)dz \\ &= \int_{[0,\infty)} p(z,t+\Delta t \,|\, x,t)\phi(z)dz \\ &\qquad + \int_{(-\infty,0]} p(z,t+\Delta t \,|\, x,t)\phi(z)dz. \end{aligned} \]

Hence, we have

\[ \int_{[0,\infty)} p_R(x,t+\Delta t \,|\, y,s) \phi(x)\,dx = \int_{[0,\infty)} \left[\int_{\mathbb{R}} p(z,t+\Delta t \,|\, x,t)\phi(z)dz \right] \,p_R(x,t \,|\, y,s)\,dx, \tag{3.23}\]

The inner integral involving \(p\) on the right-hand side of Equation 3.23 is exactly the integral we dealt with in the original derivation; we can once more Taylor expand \(\phi\) about \(z = x\), simplify and let \(\Delta t \downarrow 0\) to obtain

\[ \int_{[0,\infty)} \frac{\partial}{\partial t}p_R(x,t+\Delta t \,|\, y,s) \phi(x)\,dx = \int_{[0,\infty)} \left[\phi '(x)f(x) + \phi''(x) \frac{g(x)^2}{2} \right] p_R(x,t \,|\, y,s)\,dx. \tag{3.24}\]

We again apply integration by parts to simplify the integrals on the right-hand side. Doing so yields

\[ \begin{aligned} \int_{[0,\infty)} \phi '(x)f(x) p_R(x,t \,|\, y,s)\,dx &= - \int_{[0,\infty)}\phi(x) \frac{\partial }{\partial x}\left( f(x) p_R(x,t \,|\, y,s) \right)dx \\ &\qquad+ \phi(x)f(x) \, p_R(x,t\,|\,y,s) \vert_{x=0}, \\ \int_{[0,\infty)} \phi''(x) \frac{g(x)^2}{2} \,p_R(x,t \,|\, y,s)\,dx &= \int_{[0,\infty)} \phi(x)\frac{\partial^2 }{\partial x^2}\left( \frac{g^2(x)}{2} p_R(x,t \,|\, y,s) \right)dx \\ &\qquad -\phi(x)\left[ \frac{\partial}{\partial x}\left( \frac{g^2(x)}{2} p_R(x,t \, | \, y,s) \right) \right]_{x = 0}, \end{aligned} \]

where we used the fact that \(\phi'(0) = 0\) because \(\phi\) is an even function. Therefore, for an arbitrary test function \(\phi\), we have

\[ \begin{aligned} \int_{[0,\infty)} \frac{\partial}{\partial t}p_R(x,t+\Delta t \,|\, y,s)& \phi(x)\,dx = - \int_{[0,\infty)}\phi(x) \frac{\partial }{\partial x}\left( f(x) p_R(x,t \,|\, y,s) \right)dx \\ &+ \int_{[0,\infty)} \phi(x)\frac{\partial^2 }{\partial x^2}\left( \frac{g^2(x)}{2} p_R(x,t \,|\, y,s) \right)dx \\ &+ \left[ \phi(x)f(x) \, p_R(x,t\,|\,y,s) - \phi(x) \frac{\partial}{\partial x}\left( \frac{g^2(x)}{2} p_R(x,t \, | \, y,s) \right) \right]_{x = 0} \end{aligned} \]

The only way this equation can hold for an arbitrary \(\phi\) is if the integrands agree, i.e. \(p_R\) obeys the standard Fokker-Planck equation, and, since \(\phi(0)\) is arbitrary, we also need the last term on the right-hand side to vanish regardless of the value of \(\phi(0)\), i.e.

\[ f(x) \, p_R(x,t\,|\,y,s) - \frac{\partial}{\partial x}\left( \frac{g^2(x)}{2} p_R(x,t \, | \, y,s) \right) = 0 \quad \text{ at } x= 0. \tag{3.25}\]

Equation 3.25 is thus the condition that \(p_R\) must satisfy at \(x = 0\) in order to take into account the reflecting boundary condition on the dynamics of the underlying stochastic process \(X_R\). In fact, we can interpret condition Equation 3.25 as a no-flux boundary condition at \(x = 0\) as follows: We can write the Fokker-Planck equation Equation 3.11 in the form

\[ \frac{\partial }{\partial t} p_R(x,t \, | \, y,s) + \frac{\partial Q}{\partial x} = 0 \]

where the probability flux \(Q\) is given by

\[ Q(x) = f(x)p_R(x,t \, | \, y,s) - \frac{\partial }{\partial x} \left( \frac{g^2(x)}{2} p_R(x,t \,|\, y,s) \right). \]

The condition in Equation 3.25 corresponds to taking \(Q(0) = 0\), i.e. no flux at the boundary \(x = 0\).

A reflecting boundary is not the only possible choice of boundary condition for the Fokker-Planck equation. For example, we could decide to have an “absorbing boundary” in which we terminate trajectories that reach the boundary. This would correspond to ignoring cells which touch the edge of the Petri dish during an experiment. If we set an absorbing boundary at \(x = 0\), then the boundary condition for the Fokker-Planck equation is \(p(0,t \, | \, y,s) = 0\), i.e., there is no chance that a particle can reach \(x = 0\). This means that the integral \(\int_{(0,\infty)} p(x,t)\,dx\) is equal to \(1\) at time \(t = 0\) but its values will decay in time and obey

\[ \lim_{t \to \infty} \int_{(0,\infty)} p(x,t) \, dx = 0 \]

since all trajectories will eventually hit \(x = 0\) and be absorbed if we wait long enough.

Our choice of boundaries is typically motivated by the application we have in mind for the underlying process and other boundary conditions can be implemented using approaches similar to that employed above for the reflecting case.

3.3 The Kolmogorov Backward Equation

3.3.1 Derivation

The Fokker-Planck equation Equation 3.11 tells us how the probability density function of the process changes as time evolves, i.e., how the probability of being in a certain state changes with time. In some situations, we actually want to know how the probability of ending up in a certain state depends on where the process starts, a question that the Fokker-Planck equation is not designed to answer. Instead, for this question, we need the so-called Kolmogorov backward equation of the process. If we consider the conditional density function \(p(x,t \, | \, y,s)\), then the Fokker-Planck equation tells us what happens when we vary \(x,t\) and the Kolmogorov backward equation gives us information about what happens when \(y,s\) are varied. Moreover, we can use the backward equation to derive information about the mean hitting times of states depending on a given starting position and this will greatly improve our understanding of the noise-induced switching phenomenon we have seen in previous examples.

To derive the Kolmogorov backward equation, we begin by considering the Chapman-Kolmogorov equation with a relabelling of the variables:

\[ p(x,t\,|\, y,s - \Delta s) = \int_{\mathbb{R}} p(x,t \,|\, z,s) \,p(z,s \,|\, y,s - \Delta s)\,dz. \tag{3.26}\]

This equation is valid for any \(\Delta s>0\) and our goal will be to take the limit \(\Delta s \downarrow 0\) presently. Previously we introduced a smooth test function at this point (for regularity purposes) but this will be a more formal derivation. Thus we will take a rather risky Taylor expansion of \(p(x,t \, | \, z,s)\) about the point \(z = y\) to obtain

\[ p(x,t \, | \, z,s) = p(x,t \, | \, y,s) + (z-y) \frac{\partial }{\partial y} p(x,t \, | \, y,s) + \frac{(z-y)^2}{2} \frac{\partial^2}{\partial y^2} p(x,t \, | \, y,s) + o((z-y)^2). \]

Substituting this expression into the right-hand side of Equation 3.26 yields

\[ \begin{aligned} p(x,t\,|\, y,s - \Delta s) &= p(x,t\,|\, y,s) \int_{\mathbb{R}} p(z,s \, | \, y, s - \Delta s)\, dz \\ &\quad + \frac{\partial }{\partial y} p(x,t \, | \, y,s) \int_{\mathbb{R}} (z-y) p(z,s \, | \, y,s - \Delta s)\,dz \\ &\quad + \frac{\partial^2 }{\partial y^2} p(x,t \, | \, y,s) \int_{\mathbb{R}} \frac{(z-y)^2}{2} p(z,s \, | \, y,s - \Delta s)\,dz + o(\Delta s^2). \end{aligned} \]

The first integral on the right-hand side of this equation is \(1\) because it is the integral of a PDF. The second and third integrals are the mean and variance of the process and are dealt with as in the derivation of the Fokker-Planck equation (see equations Equation 3.17 and Equation 3.18). Simplifying these integrals thus yields

\[ \begin{aligned} p(x,t\,|\, y,s - \Delta s) &= p(x,t\,|\, y,s) + f(y)\frac{\partial }{\partial y} p(x,t \, | \, y,s) \Delta s\\ &\quad+ \frac{g(y)^2}{2}\frac{\partial^2 }{\partial y^2} p(x,t \, | \, y,s) \Delta s + o(\Delta s^2). \end{aligned} \]

Rearranging, dividing across by \(\Delta s\) and letting \(\Delta s \downarrow 0\) gives us the Kolmogorov backward equation

\[ \boxed{ -\frac{\partial }{\partial s} p(x,t \, | \, y,s) = f(y)\frac{\partial }{\partial y} p(x,t \, | \, y,s) + \frac{g(y)^2}{2}\frac{\partial^2 }{\partial y^2} p(x,t \, | \, y,s). } \] {#eq-KBE}

If we define the function \(d(y)= g^2(y)/2\), then we can write the Fokker-Planck equation more compactly as

\[ \frac{\partial }{\partial t} p(x,t) = \frac{\partial^2}{\partial x^2}\left( d(x) p(x,t) \right) - \frac{\partial }{\partial x} \left( f(x) p(x,t) \right) \]

and the Kolmogorov backward equation as

\[ -\frac{\partial }{\partial s} p(x,t \, | \, y,s) = f(y)\frac{\partial }{\partial y} p(x,t \, | \, y,s) + d(y)\frac{\partial^2 }{\partial y^2} p(x,t \, | \, y,s). \]

The stationary distribution from solving the steady-state version of the Fokker-Planck equation is then given by

\[ p_s(x) = \frac{C}{d(x)} \exp\left( \int_0^x \frac{f(y)}{d(y)} \, dy \right), \quad C > 0. \]

3.3.2 Application to mean hitting times

Consider once more the SDE from Example 3.4 earlier in this chapter:

\[ X(t + \Delta t) = X(t) + \left( -k_1 X(t)^3 + k_2 X(t)^2 - k_3 X(t) + k_4 \right)\Delta t + k_5 \sqrt{\Delta t} \, \xi. \]

Using the same parameter set as before, we know that this SDE can exhibit noise-induced switching between the stable steady states \(x_{s1} = 100\) and \(x_{s2} = 400\) (see Figure 3.2 right). There is also an unstable steady state \(x_u = 250\). Our goal will be to estimate how long on average it takes the process to hit the unstable steady state \(x_u\) given that it starts below that point. To this end, define the probability \(h(y,t)\) by

\[ h(y,t) = \mathbb{P}\left[ X(t') \in (-\infty,x_u) \,\,\forall t' < t, \,\, X(0) = y < x_u \right]. \]

In other words, \(h(y,t)\) is the probability that the process stays below \(x_u\) before time \(t\), given that it started at position \(y\). The quantity \(p(x,t \, | \, y,0)dx\) is approximately the probability that a particle is at position \(y\) at time zero and in the interval \([x,x+dx]\) at time \(t\); this density obeys the Fokker-Planck (FP) and Kolmogorov backward equations (KBEs). If we supplement the FP and KBEs with the boundary conditions

\[ p(x_u,t \, | \, y,s ) = p(x,t \, | \, x_u,s) = 0,\quad s < t, \qquad x < x_u, \]

then we ensure that \(p(x,t \, | \, y,s) = 0\) if \(y \geq x_u\) or \(x \geq x_u\). In other words, with these boundary conditions, paths of the particle that cross the level \(x_u\) before time \(t\) have probability zero. Hence we can write the desired probability as

\[ h(y,t) = \int_{-\infty}^{x_u} p(x,t \, | \, y,0)\,dx, \]

where \(p(x,t \, | \, y,0)\,dx\) is the probability that the process remains in \((-\infty,x_u)\) and lies in \([x,x+dx]\) at time \(t\), given that it started at \(y\) initially. Since \(f,g\) do not depend explicitly on time, we can then shift time as follows in the expression for \(h\):

\[ h(y,t) = \int_{-\infty}^{x_u} p(x,0 \, | \, y,-t)\,dx. \]

Similarly, shifting time in the Kolmogorov backward equation (\(s \mapsto -t\)) yields:

\[ \frac{\partial }{\partial t} p(x,0 \, | \, y,-t) = f(y)\frac{\partial }{\partial y} p(x,0 \, | \, y,-t) + d(y)\frac{\partial^2 }{\partial y^2} p(x,0 \, | \, y,-t). \]

Integrating over \(x\) thus yields

\[ \frac{\partial }{\partial t} h(y,t) = f(y)\frac{\partial }{\partial y} h(y,t) + d(y)\frac{\partial^2 }{\partial y^2} h(y,t). \tag{3.27}\]

Define the hitting time of \(x_u\) starting from \(y\) by

\[ \tau(y) = \inf\{ t > 0: \, X(t) = x_u, \, X(0)= y \}, \]

and thus the average hitting time \(\bar{\tau}(y)\) is given by \(\bar{\tau}(y) = \mathbb{E}[\tau(y)]\). Formally, we can observe that

\[ \mathbb{P}[\tau(y) = t] \approx h(y,t) - h(y,t+dt) \approx - \frac{\partial}{\partial t} h(y,t)dt. \] Next apply integration by parts to find that

\[ \bar{\tau }(y) = - \int_0^\infty t \frac{\partial }{\partial t} h(y,t) \, dt = \int_0^\infty h(y,t)\,dt. \]

Next integrate Equation 3.27 with respect to \(t\) and use the initial condition \(h(y,0) = 1\) to show that

\[ - 1 = f(y) \frac{d}{dy} \bar{\tau}(y) + d(y) \frac{d^2}{dy^2} \bar{\tau}(y), \quad y \in (-\infty,x_u). \tag{3.28}\]

Note that \(\lim_{t \to \infty} h(y,t) = 0\) because every trajectory will eventually cross \(x_u\) if we wait long enough. We need to impose some boundary conditions in order to solve Equation 3.28. If \(y = x_u\), then \(p(x,t \, | \, x_u,s) = 0\), implying that \(h(x_u,t) = 0\), which means that

\[ \bar{\tau}(x_u ) = 0. \]

Next, if we start very far from \(x_u\), then the initial condition should have a negligible impact on how long it takes to hit the level \(x_u\) and thus we assume that

\[ \frac{d}{dy}\bar{\tau}(y) \vert_{y = -\infty} = 0. \]

With these boundary conditions in hand, we can solve Equation 3.28 to show that

\[ \boxed{\bar{\tau}(y) = \int_y^{x_u} \frac{1}{d(z) p_s(z)} \int_{-\infty}^z p_s(x)\,dx\,dz.} \]

Now we can calculate, for example, how long it takes on average to hit \(x_u\) starting from the lower stable state, i.e.

\[ \bar{\tau}(x_{s1}) = \int_{x_{s1}}^{x_u} \frac{1}{d(z) p_s(z)} \int_{-\infty}^z p_s(x)\,dx\,dz. \]

Figure 3.5: Left: Paths of the SDE Equation 3.9 with different values of \(\Delta t\). Right: Estimates of the mean hitting time of the state \(x_u = 250\), \(\bar{\tau}(y)\), for different values of the initial condition \(y\). [Code: `CH3_KBE_hitting.m` plus associated .dat files]

Figure 3.5 (right) shows the results of computing the integral formula for \(\bar{\tau}\) above for different values of \(y\) and compares it to a simulation approach of estimating the hitting time over a large number of paths of the process. The left panel of Figure 3.5 demonstrates one drawback or source of inaccuracy of the direct simulation approach, namely discretisation error. The blue and red paths are essentially the same sample path but the blue path is sampled using a smaller discretisation parameter and hence should be more accurate. We observe that for this particular sample path, the red solution underestimates the hitting time of the level \(x_u = 250\). We could take smaller and smaller step sizes for accuracy but of course this will be computationally intensive and demonstrates the advantage of our analytic formula for \(\bar{\tau}(y)\).

3.4 The Chemical Fokker-Planck Equation

Thus far we have seen that SDEs can be useful modelling tools in their own right and that it is more analytically tractable to study certain problems for SDEs than it is for their chemical reaction process counterparts. For example, we developed theory to estimate mean hitting times of states for a general class of SDEs and applied this to the study of systems with noise-induced switching between alternative stable states. We will conclude this chapter by making a connection between SDEs and chemical reaction processes that will allow us to apply our analytic tools for SDEs to reaction processes, thereby enabling us to tackle some problems which were intractable up to now. Our approach will be to approximate the chemical master equations of a reaction process by an appropriate Fokker-Planck equation in the large molecule number regime; this approximation is the so-called Chemical Fokker-Planck equation.

3.4.1 Derivation for production-degradation processes

Consider our old friend the single-species production degradation process:

\[ A \xrightarrow[]{k_1} \emptyset, \quad \emptyset \xrightarrow[]{k_2} A. \tag{3.29}\]

In Chapter 1 we showed that the chemical master equations of the process Equation 3.29 are given by

\[ \frac{d}{dt}P_n(t) = k_1 (n+1)P_{n+1}(t) - k_1 n P_n(t) + k_2 \nu P_{n-1}(t) - k_2 \nu P_n(t), \quad n \geq 0, \quad t \geq 0. \]

Introducing the auxiliary functions

\[ h_1(n,t) = k_1 n P_n(t), \quad h_2(n,t) = k_2 \nu P_n(t), \]

we can rewrite the chemical master equations as

\[ \frac{d}{dt}P_n(t) = h_1(n+1,t) - h_1(n,t) + h_2(n-1,t) -h_2(n,t), \quad n \geq 0, \quad t \geq 0. \tag{3.30}\]

The key to our approximation of Equation 3.30 will be to assume that we can concentrate on the case when the number of molecules of \(A\), denoted by \(n\), is large; justifying this assumption would largely depend on the specific dynamics of the process in question. Consider a fixed large number \(\omega\) such that \(n = \eta \omega\) for some continuous variable \(\eta > 0\) and set \(P_n(t) = P(\eta,t)\). Our auxiliary functions thus become

\[ h_1(\eta,t) = k_1 \eta \omega P(\eta,t), \quad h_2(\eta,t) = k_2 \nu P(\eta,t) \]

and Equation 3.30 now reads

\[ \frac{\partial}{\partial t}P(\eta, t) = h_1\left(\eta + \tfrac{1}{\omega},t \right) - h_1(\eta,t) + h_2\left(\eta - \tfrac{1}{\omega},t \right) -h_2(\eta,t), \quad t \geq 0. \tag{3.31}\]

Given that \(\omega\) is a large parameter, we can Taylor expand the right-hand side of Equation 3.31 in powers of \(\omega\) about the point \((\eta,t)\). The expansions of \(h_1\) and \(h_2\) will be:

\[ \begin{aligned} h_1\left(\eta + \tfrac{1}{\omega},t \right) &\approx h_1(\eta, t) + \frac{1}{\omega} \frac{\partial }{\partial \eta }h_1(\eta, t) + \frac{1}{2\omega^2} \frac{\partial^2}{\partial \eta^2} h_1(\eta,t) + O(h_1/\omega^3),\\ h_2\left(\eta - \tfrac{1}{\omega},t \right) &\approx h_2(\eta, t) - \frac{1}{\omega} \frac{\partial }{\partial \eta }h_2(\eta, t) + \frac{1}{2\omega^2} \frac{\partial^2}{\partial \eta^2} h_2(\eta,t) + O(h_2/\omega^3) . \end{aligned} \]

Substituting these expansions into Equation 3.31 and truncating at \(O(1/\omega^2)\) yields:

\[ \frac{\partial}{\partial t}P(\eta, t) = \frac{1}{\omega} \frac{\partial }{\partial \eta }h_1(\eta, t) + \frac{1}{2\omega^2} \frac{\partial^2}{\partial \eta^2} h_1(\eta,t) - \frac{1}{\omega} \frac{\partial }{\partial \eta }h_2(\eta, t) + \frac{1}{2\omega^2} \frac{\partial^2}{\partial \eta^2} h_2(\eta,t). \]

For brevity and notational convenience, we can let \(x = \eta \omega\) to obtain the slightly more pleasant expression

\[ \frac{\partial}{\partial t}P(x, t) = \frac{\partial }{\partial x }h_1(x, t) + \frac{1}{2} \frac{\partial^2}{\partial x^2} h_1(x,t) - \frac{\partial }{\partial x }h_2(x, t) + \frac{1}{2} \frac{\partial^2}{\partial x^2} h_2(x,t), \tag{3.32}\]

while bearing in mind that this approximation is now only valid for large values of \(x\) by dint of the relationship between \(n\) (molecule number) and \(x\). Carefully tracing through our latest change of variables, we have \[ h_1(x,t) = k_1 x P(x,t), \quad h_2(x,t) = k_2 \nu P(x,t), \] and hence Equation 3.32 can be rewritten as

\[ \boxed{ \frac{\partial}{\partial t}P(x, t) = -\frac{\partial }{\partial x }\left( f(x) P(x,t) \right) + \frac{\partial^2}{\partial x^2}\left(d (x) P(x,t)\right), } \tag{3.33}\]

where, in our usual SDE/Fokker-Planck notation, we have

\[ \boxed{ f(x) = - k_1 x + k_2 \nu, \quad d(x) = \frac{g^2(x)}{2} = \frac{k_1 x + k_2 \nu }{2}. } \tag{3.34}\]

Thus Equation 3.33, together with the coefficients from Equation 3.34, are the chemical Fokker-Planck equation for the process Equation 3.29. The solution \(P(x,t)\) of Equation 3.33 gives us the approximate probability that we will have roughly \(x\) molecules of \(A\) present at time \(t\).

At least for large values of \(x\) (i.e., large numbers of molecules of \(A\)), we can use Equation 3.33 to answer questions about the process Equation 3.29. For example, we can solve the steady-state version of Equation 3.33 to readily obtain an approximate stationary distribution for the production-degradation process. Doing so yields

\[ \begin{aligned} P_s(x) &= \frac{C}{d(x)} \exp\left( \int_0^x \frac{f(y)}{d(y)} \, dy \right) \\ &= \frac{2 C}{k_1 x + k_2 \nu }\exp\left( 2 \int_0^x \frac{-k_1 y + k_2 \nu }{k_1 y + k_2 \nu }\,dy \right) \\ &= \frac{2 C}{k_1 x + k_2 \nu } \exp\left( -2x + 4k_2 \nu \int_0^x \frac{dy}{k_1 y + k_2 \nu} \right) \\ &= 2C \exp\left( -2x + \frac{4k_2 \nu }{k_1} \log(k_1 x + k_2 \nu) - \frac{4k_2 \nu }{k_1}\log(k_2 \nu ) \right), \quad x > 0, \quad C>0. \end{aligned} \tag{3.35}\]

The approximate stationary distribution Equation 3.35 is only valid for \(x>0\) since \(x\) represents molecule numbers in this formulation and thus the normalisation constant will be given by

\[ C = \frac{1}{2} \left[ \int_0^\infty \exp\left( -2x + \frac{4k_2 \nu }{k_1}\log(k_1 x + k_2 \nu) - \frac{4k_2 \nu }{k_1}\log(k_2 \nu ) \right)\,dx \right]^{-1}. \]

Figure 3.6 (left) shows excellent agreement between the stationary distribution Equation 3.35 from the chemical Fokker-Planck approximation and estimated stationary distribution obtained via direct simulation of the process Equation 3.29.

Figure 3.6: Left: Comparison of the stationary distribution Equation 3.35 and estimated stationary distribution from simulations of the process Equation 3.29. Right: Average hitting time of 19 molecules as a function of the initial condition estimated from direct simulations and the integral formula. Parameters: \(k_1 = 0.1\), \(k_2\nu = 1\). [Code(s): `CH3_CFP_HT1.m` and `CH3_CFP_HT2.m`]

Armed with an (approximate) analytic expression for the stationary distribution, we can then turn to estimating mean hitting times for the process Equation 3.29. For example, we can define

\[ \tau(y) = \inf\{ t > 0: \, A(t) = 19, \,\, A(0) = y \}, \quad y \in \{0,1,\dots,18\}. \]

From the theory developed earlier in this chapter, we have

\[ \mathbb{E}[\tau(y)] = \int_y^{19} \frac{1}{d(z) P_s(z)}\int_0^z P_s(x)\,dx\,dz, \] where \(P_s\) is given by Equation 3.35 and the function \(d\) is given by Equation 3.34. Figure 3.6 (right) shows \(\mathbb{E}[\tau(y)]\) as a function of the initial condition (\(y\)) compared to the average hitting time obtained via direct simulation of the process. Once more, theory and numerical experiment are in very close agreement with the chemical Fokker-Planck equation offering an accurate approximation without the computational overhead of direct simulations.

3.4.2 General multi-species framework

Before we can write the general chemical Fokker-Planck equation we need to write down the multi-dimensional analogue of the Fokker-Planck equation Equation 3.11. To this end, consider the \(N\)-dimensional SDE

\[ \textbf{X}(t+\Delta t) = \textbf{X}(t) + \textbf{f}(\textbf{X}(t))\Delta t + \textbf{g}(\textbf{X}(t)) \sqrt{\Delta t}\, \boldsymbol{\xi}, \quad t > 0, \tag{3.36}\]

where \(\textbf{X}(t)\) and \(\textbf{f}(\textbf{X}(t))\) are \(N\)-dimensional vectors, \(\textbf{g}(\textbf{X}(t))\) is an \(N \times N\) matrix and \(\boldsymbol{\xi}\) is an \(N\)-dimensional vector of independent standard Normal random variables. We define the diffusion tensor to be the matrix \(\textbf{D} = \tfrac{1}{2} \textbf{g}(x) \textbf{g}^T(x)\), i.e.,

\[ D_{i,j}(\textbf{x}) = \frac{1}{2}\sum_{k = 1}^N g_{i,k}(\textbf{x})g_{j,k}(\textbf{x}). \]

The Fokker-Planck equation for the \(N\)-dimensional process \(\textbf{X}\) solving Equation 3.36 is given by

\[ \frac{\partial }{\partial t}P(\textbf{x},t) = - \sum_{i = 1}^N \frac{\partial }{\partial x_i} \left[ f_i(\textbf{x})P(\textbf{x},t) \right] + \sum_{i = 1}^N \sum_{j=1}^N \frac{\partial^2 }{\partial x_i \partial x_j} \left[ D_{i,j}(\textbf{x}) P(\textbf{x},t) \right]. \]

Now suppose that we have a chemical reaction process in which:

we have \(N \geq 1\) chemical species that can undergo \(q \geq 1\) distinct reactions,
the \(N\)-dimensional vector \(\textbf{X}(t)\) will represent the state of the system at time \(t\), i.e., \(X_i(t)\) is the number of molecules of species \(i\) at time \(t\),
\(\alpha_j(\textbf{x})\) is the propensity function associated with the \(j\)th reaction, where \(\textbf{x} \in \mathbb{R}^N\) is the current state of the system,
\(v_{j,i}\) is the change in the number of molecules in species \(i\), i.e., \(X_i\), that occurs when reaction \(j\) takes place.

The chemical master equations for this reaction process can be written in the form:

\[ \frac{\partial }{\partial t} P(\textbf{x},t) = \sum_{j = 1}^q\left( \alpha_j(\textbf{x} - \textbf{v}_j)P(\textbf{x} - \textbf{v}_j, t) - \alpha_j(\textbf{x})P(\textbf{x},t) \right), \]

where \(P(\textbf{x},t) = \mathbb{P}[\textbf{X}(t) = \textbf{x}]\).

Carrying out the same large molecule number approximation as we did for the production-degradation example, we can obtain the multi-dimensional chemical Fokker-Planck equations:

\[ \begin{aligned} \frac{\partial }{\partial t} P(\textbf{x},t) &= - \sum_{i = 1}^N \sum_{j=1}^q \frac{\partial }{\partial x_i} \left[ v_{j,i}\alpha_j(\textbf{x})P(\textbf{x},t) \right] \\ &\qquad+ \frac{1}{2}\sum_{i = 1}^N \sum_{k=1}^N \frac{\partial^2 }{\partial x_i \partial x_k} \left[\sum_{j=1}^q v_{j,i}v_{j,k} \alpha_j(\textbf{x}) P(\textbf{x},t)\right]. \end{aligned} \]

Knowledge checklist

Key topics:

Stochastic differential equations (computational definition and SSAs)
The Fokker-Planck equation (stationary distribution, boundary conditions)
The Kolmogorov Backward equation (expected hitting times)
The Chemical Fokker-Planck equation (approximating chemical reaction processes)

Key skills:

Use the computational definition of an SDE to compute moments and other related quantities.
State and derive the Fokker-Planck equation for a given SDE, including applying the appropriate boundary conditions.
Calculate stationary distributions of SDEs using the steady-state Fokker-Planck equation.
Compute hitting times for SDEs.
Explain the impact of boundary conditions on stationary distributions and hitting times.
Explain how SSAs for SDEs can be modified to account for a range of boundary conditions that occur in biological applications (e.g., reflecting boundaries).
Write chemical Fokker-Planck approximation for a given chemical reaction process.